######## INFO ########

# PROJECT
# Replication files for: Economic Impact of the South China Sea Dispute in China-Philippines
#                        Relations from 2012 to 2016: A DM-LFM analysis  

# Authors: Dianyi Yang
# R Script
# Purpose: This script produces the remaining figures for imports from China, 
#          and remaining tables.
# Created: 30 Oct 2023
# Updated: 29 Nov 2023
# Inputs: stored/dmlfm_import.rdata
#         custom_functions/DM-LFM.r
#         stored/dmlfm.rdata
#         stored/dmlfm_placebo_import.rdata
#         stored/dmlfm_placebo.rdata
#         stored/leave-one-out.rdata
# Outputs: output/dmlfm_trend_import.pdf
#          output/dmlfm_trend_import.png
#          output/dmlfm_effect_import.pdf
#          output/dmlfm_effect_import.png
#          Table 4: Main Analysis Result table for the DM-LFM.
#          Table 5: Result table for the ATT over time
#          output/placebo_trend_import.pdf
#          output/placebo_trend_import.png
#          Table 6: Placebo result table for the DM-LFM.
#          output/loo_trend_import.pdf
#          output/loo_trend_import.png
#          output/loo_effect_import.pdf
#          output/loo_effect_import.png
#          Table 7: Result table for the ATT estimated by SCM.

######## SETUP ########
rm(list = ls()) # clear workspace

need <- c("Synth", "reshape2","devtools", "synthdid",'tidyverse',
          'bpCausal', "fixest","modelsummary",
          'kableExtra') # list packages needed
have <- need %in% rownames(installed.packages()) # checks packages you have
if(any(!have)) install.packages(need[!have]) # install missing packages
#devtools::install_github("synth-inference/synthdid") #manually install sdid if the codes above fail to do so automatically
invisible(lapply(need, library, character.only=T)) # load needed packages

setwd(gsub('/code','',dirname(rstudioapi::getSourceEditorContext()$path))) #set wd to root directory
list.files()
##################DM-LFM####################################
load('stored/dmlfm_import.rdata')
source('custom_functions/DM-LFM.r')

result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
#plot.DMLFM(out, ci.level=0.90)
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
class(result.use90) <- 'DMLFM'

#trend plot
ggplot(result.use$est.eff, aes(x=time))+
  geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "95% CI"), alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "90% CI"), alpha = 0.5)+
  geom_line(aes(y=observed, color='Observed'),linewidth=0.75,alpha=0.5)+
  geom_line(aes(y=estimated_counterfactual, color='Counterfactual'),linewidth=0.75,alpha=0.5)+
  geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+ 
  scale_color_manual(name= '',
                     breaks=c('Observed','Counterfactual'),
                     values=c('Observed'="#00BFC4", "Counterfactual"="#F8766D"))+
  labs(x="Months since treatment", y='Log(Filipino imports from China)')+
  scale_fill_manual(name=c("", ''),
                    breaks=c('95% CI', '90% CI'),
                    values=c('lightgrey', 'grey'))  
ggsave('output/dmlfm_trend_import.pdf', width = 10, height=5) #Figure 4
ggsave('output/dmlfm_trend_import.png', width = 10, height=5)

#effect plot
ggplot(result.use$est.eff, aes(x=time, y=estimated_ATT))+
  geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u, fill = "95% CI"), alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u, fill = "90% CI"), alpha = 0.5)+
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
  theme_bw()+labs(x="Months since treatment",y="ATT with 90% and 95% CIs") + geom_line(linewidth=0.75) +
  scale_fill_manual(name=c("", ''),
                    breaks=c('95% CI', '90% CI'),
                    values=c('lightgrey', 'grey'))         #effect plot
ggsave('output/dmlfm_effect_import.pdf', width = 10, height=5) #Figure 6
ggsave('output/dmlfm_effect_import.png', width = 10, height=5)

#Regression Table
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'

result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'

result.use_import <- result.use
result.use90_import <- result.use90

load('stored/dmlfm.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI

result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'

result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'


tidy.DMLFM <- function(x, ...) {
  ret <- data.frame(
    term      = 'ATT',
    estimate  = x$est.avg[[1]],
    conf.low  = x$est.avg[[2]],
    conf.high = x$est.avg[[3]],
    p.value   = x$est.avg_p)
  ret
}

glance.DMLFM <- function(x, ...) {
  ret <- data.frame(
    nobs  = x$obs,
    tr  = x$tr,
    ct  = x$ct)
  ret
}

gof_map <- list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0),
                list("raw" = "tr", "clean" = "Treated Units", "fmt" = 0),
                list("raw" = "ct", "clean" = "Control Units", "fmt" = 0)) 
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90, #Table 4
                  '95% CI'=result.use_import, '90% CI'=result.use90_import),
             estimate = "{estimate}", #output='latex', #remove hashtag for latex output
             statistic = "[{conf.low}, {conf.high}]", 
             gof_map = gof_map, escape=T,
             notes=c('Equal-tailed Credible Intervals in square brackets.'),
             title = 'Main Analysis Result table for the DM-LFM.')  %>% 
  add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)


#Short term result table
class(result.use) <- 'DMLFMtime'
class(result.use90) <- 'DMLFMtime'
class(result.use_import) <- 'DMLFMtime'
class(result.use90_import) <- 'DMLFMtime'
tidy.DMLFMtime <- function(x, ...) {
  ret <- data.frame(
    term      = x$est.eff$time,
    estimate  = x$est.eff$estimated_ATT,
    conf.low  = x$est.eff$estimated_ATT_ci_l,
    conf.high = x$est.eff$estimated_ATT_ci_u)
  ret
}

glance.DMLFMtime <- function(x, ...) {
  ret <- data.frame(
    nobs  = x$obs,
    tr  = x$tr,
    ct  = x$ct)
  ret
}
gof_map <- list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0),
                list("raw" = "tr", "clean" = "Treated Units", "fmt" = 0),
                list("raw" = "ct", "clean" = "Control Units", "fmt" = 0))
coef_map <- list('0' = '0',
                 '1' = '1',
                 '2' = '2',
                 '3' = '3',
                 '4' = '4',
                 '5' = '5',
                 '6' = '6')
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90, #Table 5
                  '95% CI'=result.use_import, '90% CI'=result.use90_import), 
             estimate = "{estimate}", #output='latex', #remove hashtag for latex output
             statistic = "[{conf.low}, {conf.high}]",
             gof_map = gof_map, coef_map=coef_map, escape=T,
             notes=c('Equal-tailed Credible Intervals in square brackets.'),
             title = 'Result table for the ATT over time.') %>% 
  add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)


#Placebo test (DMLFM)
load('stored/dmlfm_import_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI

#trend plot
ggplot(result.use$est.eff, aes(x=time-25))+
  geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "95% CI"), alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "90% CI"), alpha = 0.5)+
  geom_line(aes(y=observed, color='Observed'),linewidth=0.75,alpha=0.5)+
  geom_line(aes(y=estimated_counterfactual, color='Counterfactual'),linewidth=0.75,alpha=0.5)+
  geom_vline(xintercept=-25)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+ 
  scale_color_manual(name= '',
                     breaks=c('Observed','Counterfactual'),
                     values=c('Observed'="#00BFC4", "Counterfactual"="#F8766D"))+
  labs(x="Months since actual treatment", y='Log(Filipino imports from China)')+
  scale_fill_manual(name=c("", ''),
                    breaks=c('95% CI', '90% CI'),
                    values=c('lightgrey', 'grey'))  
ggsave('output/placebo_trend_import.pdf', width = 10, height=5) #Figure 8
ggsave('output/placebo_trend_import.png', width = 10, height=5)

#effect plot (90% and 95% CIs)
ggplot(result.use$est.eff, aes(x=time-25, y=estimated_ATT))+
  geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), fill = "lightgrey", alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), fill = "grey", alpha = 0.5)+
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  geom_vline(xintercept = -25, linetype = "dashed", color = "grey") +
  geom_line(linewidth=0.75)+
  theme_bw()+labs(x="Months since treatment",y="ATT with 90% and 95% CIs")

#placebo table
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'

result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'

result.use_import <- result.use
result.use90_import <- result.use90

load('stored/dmlfm_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI

result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'

result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'


modelsummary(list('95% CI'=result.use, '90% CI'=result.use90, #Table 6
                  '95% CI'=result.use_import, '90% CI'=result.use90_import), 
             estimate = "{estimate}", #output='latex', #remove hashtag for latex output
             statistic = "[{conf.low}, {conf.high}]", 
             gof_map = gof_map, escape=T,
             notes=c('Equal-tailed Credible Intervals in square brackets.'),
             title = 'Placebo result table for the DM-LFM.') %>%
  add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)

#leave-one-out test
load('stored/leave-one-out_import.rdata')

narrow <- plot %>%
  group_by(time) %>%
  summarize(counterfactual_ci_l = max(counterfactual_ci_l), counterfactual_ci_u = min(counterfactual_ci_u))

#trend
ggplot(plot, aes(x=time))+
  geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u), alpha = 0.5, fill='lightgrey')+
  geom_ribbon(data=narrow,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u), alpha = 0.5, fill='grey')+
  geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_counterfactual, color= miss_unit),linewidth=0.5,alpha=0.5)+
  geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_counterfactual, color='Counterfactual (full pool)'),linewidth=1,alpha=1)+
  geom_line(aes(y=observed, color='Observed'),linewidth=1,alpha=1)+
  geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
  labs(x="Months since treatment", y='Log(Filipino imports from China)')+
  scale_color_discrete(name='Key',
                       breaks=c('Observed', 'Counterfactual (full pool)', unique(plot$miss_unit)))
ggsave('output/loo_trend_import.pdf', width = 10, height=5) #Figure 10
ggsave('output/loo_trend_import.png', width = 10, height=5)

#effect
narrow <- plot %>%
  group_by(time) %>%
  summarize(estimated_ATT_ci_l = max(estimated_ATT_ci_l), estimated_ATT_ci_u = min(estimated_ATT_ci_u))

ggplot(plot, aes(x=time))+
  geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='lightgrey')+
  geom_ribbon(data=narrow,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='grey')+
  geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_ATT, color= miss_unit),linewidth=0.5,alpha=0.5)+
  geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_ATT, color='ATT with full pool'),linewidth=1,alpha=1)+
  geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
  labs(x="Months since treatment", y='ATT with widest and narrowest 95% CIs')+
  scale_color_discrete(name='Key',
                       breaks=c('ATT with full pool', unique(plot$miss_unit)))
ggsave('output/loo_effect_import.pdf', width = 10, height=5) #Figure 12
ggsave('output/loo_effect_import.png', width = 10, height=5)

#SC result table
load('stored/sc_import.rdata')
sc_estimate_import <- sc_estimate
load('stored/sc.rdata')
tidy.sc_estimate <- function(x, ...) {
  ret <- data.frame(
    term      = 'ATT',
    estimate  = x$estimate,
    std.err   = x$se,
    p.value   = x$p,
    conf.low  = x$estimate+x$se*qnorm(0.025),
    conf.high  = x$estimate+x$se*qnorm(0.975))
  ret
}

glance.sc_estimate <- function(x, ...) {
  ret <- data.frame(
    nobs  = 891,
    tr  = 1,
    ct  = 8)
  ret
}
gof_map <- list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0),
                list("raw" = "tr", "clean" = "Treated Units", "fmt" = 0),
                list("raw" = "ct", "clean" = "Control Units", "fmt" = 0))
modelsummary(list('Exports to China' = sc_estimate, 'Imports from China' = sc_estimate_import), 
             estimate = "{estimate}{stars} ({p.value})",
             statistic = "[{conf.low}, {conf.high}]", #output='latex', #remove hashtag for latex output
             gof_map = gof_map, escape=T,
             notes=c('p-value in parentheses; 95% CI in square brackets.'),
             title = 'Result table for the ATT estimated by SCM.') #Table 7